---
title: "Sustainable food systems"
author: "Harold Achicanoy"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    social: menu
    source_code: embed
---

```{r setup, include=FALSE}
suppressMessages(library(highcharter))
suppressMessages(library(dplyr))
suppressMessages(library(viridisLite))
suppressMessages(library(forecast))
suppressMessages(library(treemap))
suppressMessages(library(flexdashboard))
suppressMessages(library(VIM))
suppressMessages(library(plspm))
suppressMessages(library(FactoMineR))
suppressMessages(library(networkD3))
suppressMessages(library(igraph))
suppressMessages(library(plyr))
suppressMessages(library(cluster))
suppressMessages(library(plotly))

thm <- 
  hc_theme(
    colors = c("#1a6ecc", "#434348", "#90ed7d"),
    chart = list(
      backgroundColor = "transparent",
      style = list(fontFamily = "Source Sans Pro")
    ),
    xAxis = list(
      gridLineWidth = 1
    )
  )

```

Column {data-width=600}
-----------------------------------------------------------------------

### Sales Forecast

```{r}
AirPassengers %>% 
  forecast(level = 90) %>% 
  hchart() %>% 
  hc_add_theme(thm)
```

### Sustainable food system index by country

```{r}
data("worldgeojson")

indices <- readRDS(file = "//dapadfs/Workspace_cluster_9/Sustainable_Food_System/Results/sfs_index_knn_imputed.RDS")
indices[,1:(ncol(indices)-1)] <- round(indices[,1:(ncol(indices)-1)], 2)

n <- 4
colstops <- data.frame(
  q = 0:n/n,
  c = substring(viridis(n + 1), 0, 7)) %>%
  list.parse2()

highchart(type = "map") %>%
  hc_add_series_map(map = worldgeojson, df = indices, value = "SUFS", joinBy = "iso3") %>%
  hc_colorAxis(stops = color_stops()) %>%
  hc_tooltip(useHTML = TRUE, headerFormat = "",
             pointFormat = "{point.name} has a SFS index of {point.SUFS}") %>%
  hc_colorAxis(stops = colstops) %>%
  hc_legend(valueDecimals = 0, valueSuffix = "%") %>%
  hc_mapNavigation(enabled = TRUE) %>%
  hc_add_theme(thm)
```

Column {.tabset data-width=400}
-----------------------------------------------------------------------

### PLS-PM 1 - Estimations

```{r}
complete_data <- readRDS(file = "//dapadfs/Workspace_cluster_9/Sustainable_Food_System/Input_data/data_joined.RDS")

# Method 1: k nearest neighbors (non-parametric alternative)
complete_data1 <- VIM::kNN(data = complete_data); complete_data1 <- complete_data1[,colnames(complete_data)]

# PLS-PM: Using repeated indicators
# Define path model matrix (inner model)
NUTR <- c(0, 0, 0, 0)
HINT <- c(0, 0, 0, 0)
FSCY <- c(0, 0, 0, 0)
SUFS <- c(1, 1, 1, 0)
sfs_path <- rbind(NUTR, HINT, FSCY, SUFS); rm(NUTR, HINT, FSCY, SUFS)
colnames(sfs_path) <- rownames(sfs_path)
sfs_blocks <- list(2:3, 4:5, 6:9, 2:9)
sfs_modes <- rep("A", 4)
sfs_pls <- plspm(complete_data1, sfs_path, sfs_blocks, modes = sfs_modes)
plot(sfs_pls)
```

### PLS-PM 1 - Indices correlations

```{r}
pairs(sfs_pls$scores, pch = 20)
```

### PLS-PM 1 - Crossloadings analysis

```{r}
xloads = melt(sfs_pls$crossloadings, id.vars = c("name", "block"))
gg <- ggplot(data = xloads, aes(x = name, y = value, fill = block))
gg <- gg + geom_hline(yintercept = 0, color = "gray75")
gg <- gg + geom_hline(yintercept = c(-0.5, 0.5), color = "gray70", linetype = 2)
gg <- gg + geom_bar(stat = 'identity', position = 'dodge')
gg <- gg + facet_wrap(block ~ variable)
gg <- gg + theme(axis.text.x = element_text(angle = 90), line = element_blank())
ggplotly(gg)
```

### Network analysis

```{r}
rownames(complete_data1) <- complete_data1$ISO3

# Calculate similarity measure
sfs_dis <- cluster::daisy(x = complete_data1[,-1], metric = c("gower"), stand = FALSE)
sfs_dis <- 1 - as.matrix(sfs_dis)

# Do cluster analysis
sfs_pca <- FactoMineR::PCA(X = complete_data1[,-1], scale.unit = T, graph = F)
sfs_hpc <- FactoMineR::HCPC(res = sfs_pca, nb.clust = -1, graph = F)
complete_data1$cluster <- sfs_hpc$data.clust$clust

# Visualize using networkD3
sfs_dis[lower.tri(sfs_dis, diag = TRUE)] <- NA
sfs_dis <- na.omit(data.frame(as.table(sfs_dis))); names(sfs_dis) <- c("from", "to", "similarity")
sfs_dis <- sfs_dis[sfs_dis$similarity >= .98,] # Filter by more than 98 degree of similarity

gD <- igraph::simplify(igraph::graph.data.frame(sfs_dis, directed = FALSE))
nodeList <- data.frame(id = c(0:(igraph::vcount(gD) - 1)), name = igraph::V(gD)$name) # because networkD3 library requires IDs to start at 0
getNodeID <- function(x){ which(x == igraph::V(gD)$name) - 1 } # to ensure that IDs start at 0
edgeList <- plyr::ddply(sfs_dis, .variables = c("from", "to", "similarity"), 
                        function (x) data.frame(fromID = getNodeID(x$from), 
                                                toID = getNodeID(x$to)))
nodeList <- cbind(nodeList, nodeDegree = igraph::degree(gD, v = igraph::V(gD), mode = "all")); rm(gD, getNodeID)
nodeList$cluster <- as.numeric(as.character(complete_data1$cluster))[match(nodeList$name, complete_data1$ISO3)]

networkD3::forceNetwork(Links = edgeList,
                        Nodes = nodeList,
                        Source = "fromID",
                        Target = "toID",
                        Value = "similarity",
                        NodeID = "name",
                        Group = "cluster",
                        opacity = 1,
                        fontSize = 15)
```